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The evolution of dark matter in central areas of galaxies is considered (the Milky Way is taken as an example). 
It is driven by scattering off of dark matter particles by bulge stars, their absorption by the supermassive black 
hole and self-annihilation. This process is described by diffusion equation in the phase space of energy and 
angular momentum. The equation was integrated for several different models of initial dark matter distribution 
and using various assumptions about the dynamical factors. It turns out that because the Milky Way center 
is rather dynamically old (~ 4 relaxation times t r ), the difference in initial conditions almost vanishes. The 
density attains a nearly universal profile, and the 7-ray flux from dark matter annihilation lies in rather narrow 
range, which enables more robust determination of the dark matter parameters. By present time the mass of 
dark matter inside the black hole sphere of influence (r < 2 pc) has been reduced approximately twice, mostly 
because of heating by stars. It is shown that the dynamics of dark matter for t > t r is determined mainly by 
stars outside the sphere of influence. 
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I. INTRODUCTION 

It is widely accepted that the largest fraction of matter in 
the Universe is the dark matter (DM) |Q]]. Its nature is still 
unclear, but most likely it consists of yet undiscovered cold 
(non-relativistic) particles, which have very low interaction 
cross-sections with each other and with baryonic matter J2]. 
So the main physical mechanism governing the evolution of 
dark matter is gravitation. 

The dark matter is responsible for formation of large-scale 
structure of the Universe, as well as extended DM haloes of 
galaxies yfl. However, the very centers of galaxies (except 
low surface brightness galaxies) are dominated by baryons, 
which form galactic bulges. Their influence on DM haloes 
consists of at least two effects. Firstly, as the baryons cool and 
settle down in the center of a potential well created by DM, 
they change the common gravitational potential, which leads 
to compression of DM halo. This effect is called baryonic 
compression and increases DM density several times 0, H 
@, 0]. On the other side, dark matter particles are scattered 
off by bulge stars and captured by supermassive black holes 
(SMBHs), which reduces their density in the inner parts of 
bulges QS|, |9l LlO|] . Additionally, a non-zero annihilation cross- 
section of DM particles also leads to decrease of DM density 
in the very centers of galaxies and gives possibility of indirect 
detection of DM through its annihilation radiation 

ED, EH 

14]. All these issues are addressed in the present paper. 

The paper is organized as follows. In the second section we 
describe previous investigations on this subject and explain 
motivation of this study. In the third section we present our 
class of initial DM halo models, which is somewhat broader 
than in previous papers. Then the kinetic equation governing 
the evolution of DM is derived, along with its coefficients and 
boundary conditions. The fourth section is devoted to solu- 
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tion of the above equation, given for a series of initial halo 
models with inclusion of different physical processes and ap- 
proximations. This helps us to clarify validity of some com- 
monly used assumptions, and determine parameters on which 
the evolution depends mostly. Predictions for annihilation ra- 
diation are also given in this section. Finally, the conclusions 
are presented. 

II. OVERVIEW OF DARK MATTER EVOLUTION IN 
GALACTIC CENTERS 

The question about dark matter content in central areas of 
galaxies has long been studied. First of all, it depends on ini- 
tial structure of DM halo. There are two principal classes of 
methods of halo structure modelling. Analytical methods gen- 
erally consider the nonlinear stage of collapse of initial over- 
density region. These include kinetic theory of phase mixing 
lfl5ll or variants of self-similar infall models llq, ITtI [l8tl . In 
general, they predict a power-law density profile in the center 
of a DM halo, 

p d oc r~ ld , (1) 

where the power-law index 7^ is between 1 and 2. A sig- 
nificant drawback of this approach is very limited ability to 
handle mergers, especially major ones. Alternatively, numer- 
ical simulations of dissipationless A^-body gravitational inter- 
action have been largely used in recent years. They also gen- 
erally predict density profiles that have power-law cusps in 
central areas (e.g. Navarro, Frenk, White (NFW) profile lfl9tl 
with 7d = 1 or Moore profile I20I1 with 7^ = 1.5). Other 
models were proposed in which the density slope tends to zero 
at r — > (Einasto profiles 1I2TI l22tl or Burkert profiles l23ll ). 
There is still no consensus on the question whether dark mat- 
ter halos are cuspy or cored in the centers, which means that 
central density profiles obtained in simulations may equally 
good be described with different fitting formulae, although 
the inferred densities in the central parsecs would vary orders 
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of magnitude. The largest simulations of a galaxy-sized halo 
(e.g. II24|0 provide measurements of density no closer to the 
galaxy center than 3 x 10 — 3 rw — 10 3 P c ^ which is still many 
orders of magnitude greater than the distances relevant to an- 
nihilation and scattering. 

However, there are other processes that modify dark matter 
distribution in the galactic centers after halo formation. First 
of all, it is now well established that virtually all galaxies har- 
bour supermassive black holes (SMBHs) in their centers l25ll . 
whose masses Mbh range from 10 6 M Q to more than 10 9 M Q . 
A black hole (BH) immersed in a stellar cluster (or bulge as 
a whole) modifies gravitational potential at distances smaller 
than the gravitational influence radius of the black hole Th, 
defined as containing stellar mass equal to 2Mbh- If the black 
hole grew adiabatically, the surrounding matter (stars or DM) 
forms a "spike" with power-law index 7' = (9 — 2"f)/(4 — 7), 
7 being the power-law index of density profile before black 
hole formation id. 

The above argument applies for spherically symmetric evo- 
lution. If the black hole formed off-center and then spiralled 
in, the resulting profile should be shallower lf26ll . Addition- 
ally, in the hierarchical merging scenario, binary SMBHs 
should be ubiquitous. A binary SMBH ejects matter from 
r < Th due to slingshot effect, and after coalescence the stellar 
and dark matter density around it is reduced [27]. However, it 
is unlikely that a major merger took place in the recent 10 10 yr 
of Milky Way history IToll . 

If the DM consists of self-annihilating particles, the product 
of annihilation cross-section a a and relative velocity v r being 
weakly dependent on v r , then the maximum density after a pe- 
riod of time t is believed to be set at the so-called annihilation 
plateau: 

Pa=m p /(a a v r )t (2) 

(m p is the particle mass) However, this result is exact only 
for circular particle orbits. If we consider arbitrary velocity 
anisotropy ft = 1 — o\ /2<7 2 (particulary, isotropic velocity 
corresponding to ft = 0), then we obtain a weak cusp instead 
of a plateau: its density varies as r~C +1 ' 2 ' and approximately 
equals corresponding plateau density p a at the radius where 
plateau itself should transform into initial density profile. This 
issue is explained in detail in 12811 . 

Finally, a very important process that modifies the dark 
matter distribution is the gravitational scattering of DM par- 
ticles by stars in the vicinity of the SMBH and capture of par- 
ticles by the SMBH. The detailed analysis of this process is 
the main subject of the present investigation. 

The stellar distribution around a massive black hole was 
much studied in 1970s ll29l l30l l3lll . The usual method in- 
cludes representation of the distribution function (DF) / in 
terms of energy E and angular momentum L per unit mass 
(for spherically symmetric problem these are the only two 
independent variables) and consideration of diffusion in the 
{E, L} phase space. In these early studies it was shown that 
diffusion along L axis leads to formation of logarithmic de- 
pendence of / on L, so that there is continuous flux of stars to- 
wards low values of L, where they are captured by black hole. 
The diffusion along E axis drives the DF to a quasi-stationary 



state / oc |i?| 1//4 , when the loss of stars is balanced by a in- 
ward flux from higher values of E (outside the BH sphere 
of influence). The corresponding density profile p oc r~ 7 / 4 is 
called Bahcall-Wolf cusp, and is indeed observed in the center 
of our Galaxy ll32ll33i[34ll . though the density slope is some- 
what lower. It is important that the stationary solution applies 
only to systems older than two-body relaxation time ll35ll 

t _ 0.34cr 3 

G 2 m*p it In A ' 

where a is 1-d velocity dispersion, p+ - stellar density, In A w 
15 - Coulomb logarithm. It appears that in the center of our 
Galaxy t r ~ 2.5 x 10 9 yr < t Hubble and weakly depends on 
radius at r < r/j. (We adopt Mbh — 3 x 10 6 M© for consis- 
tency with other papers, though recent estimates are somewhat 
higher [36]. This yields w 2 pc). 

Dark matter dynamics is governed by the same mechanism 
of relaxation on stars, but due to negligible particle mass the 
effect of dynamic friction is unimportant, and resulting den- 
sity profile should be pd oc r~ 3 / 2 . This simple argument 
was considered in |[37ll . A more elaborate treatment involves 
time-dependent solution of diffusion equation. This was first 
done in ]8|, |9J] for a constrained initial DM distribution func- 
tion using only diffusion along angular momentum. In j38tl 
the problem was treated as one-dimensional diffusion for en- 
ergy with inclusion of loss-term due to angular momentum 
diffusion and capture by BH, which was taken from ll3lll . In 
subsequent papers llciloi the effect of self-annihilations was 
also included. These studies have led to the conclusion that 
the scattering off by stars drives DM distribution towards a 
quasi-stationary state with diminishing amplitude of density 
profile due to combined effect of heating by stars and capture 
by the SMBH. This density profile is regenerated within one 
tr even after a cusp has been destroyed by a binary SMBH 
& 

We want to draw attention to several important issues that 
were not thoroughly investigated in previous studies. 

i. The relaxation time actually depends on angular mo- 
mentum L, as will be shown later. For small values 
of L the diffusion along energy goes faster, so that the 
two-dimensional problem cannot be reduced to one- 
dimensional. 

ii. The black hole mass might have grown in time, which 
moves up the radius of BH influence and changes relax- 
ation timescale inside this radius. 

iii. The initial DM distribution function need not be 
isotropic, i.e. might depend on angular momentum. 

All these issues are examined in the present paper. Addition- 
ally, the following features are discussed: 

i. A broad class of initial distribution functions is consid- 
ered. This enables direct comparison between different 
studies that took different initial conditions, and helps 
to disentangle the effect of initial conditions from other 
effects investigated. 
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ii. The diffusion is considered both inside and outside the 
SMBH sphere of influence. Values of diffusion coeffi- 
cients are different in these two domains. 

iii. Analytical approximations of evolution of DM density 
profile and distribution function are made for limiting 
cases of t <C t r and t 3> t r . They are in good agreement 
with the solution of exact equation. 

iv. The two-dimensional diffusion equation is solved di- 
rectly. The results are compared to the simplified case 
of one-dimensional treatment, to validate the possibility 
of such simplification. 

v. The effect of self-annihilations is included as well. The 
proportion of captured, evaporated and annihilated dark 
matter mass is calculated. 

Having stated these important goals, we proceed to the in- 
vestigation itself. 

HI. PROBLEM DEFINITION 

A. Initial dark matter halo models 

Plenty of models have been suggested for DM haloes. As 
we are interested in dark matter dynamics in galactic centers, 
we shall consider only internal part of halo, which can be rea- 
sonably approximated by a power-law density profile £[). We 
operate with a distribution function f(E, L) so we need to 
transform pd to f(E, L). This transformation is essentially 
non-unique - one can obtain different DFs for the same den- 
sity profile, which differ in dependence on angular momentum 
L. It is convenient to choose a factorized family of models 
f(E, L) = fi(E) f2(R), where we have changed L for a di- 
mensionless scaled quantity Re [0..1]: 

R = L 2 /L 2 C {E) , (4) 

L c being angular momentum of a circular orbit with given en- 
ergy. (This has been suggested in A3 ill ). In the case of power- 
law density fi(E) should also be power-law 113511 . 

We consider two family of DFs, which differ in dependence 
on L: 

Model A: f(E, R) = f E 1/2 5(R ~ R ) , (5a) 
Model B: f{E, R) = f £ji/2-(i-/3)(4- 7 <i)/(2-7 d ) ir* 3 .(5b) 

Model A has been proposed in lfl5ll as a result of analytical 
treatment of DM halo collapse and phase mixing. The quan- 
tity Rq is related to orbital eccentricities of particles and is 
small, Rq <C 1. (Further discussion can be found in 10]). 
Model B is the simplest generalization of isotropic distribu- 
tion for the case of arbitrary radial velocity anisotropy f3, con- 
sidered in [40]. Here (3 = 1 — a 2 /2a 2 is Binney's anisotropy 
parameter [35]: systems with (3 = are isotropic, with 
< j3 < 1 - have radially biased velocities. For realistic 
DM haloes, (3 in the center is zero or slightly positive 11711 . as 



TABLE I: Initial halo models used in calculations, yd is the slope 
of density profile in the center before baryonic compression, f3 is 
the velocity anisotropy parameter, ph and p' h are density values at 
r = r;, = 2 pc before and after compression, in Afo/pc 3 . 



Model 


Id 


13 


Ph 


Ph 


Al 


1.7 


0.6 


10 4 


2.3 x 10 4 


Bl 


1.5 





2 x 10 3 


1.3 x 10 4 


B2 


1.5 


0.5 


2 x 10 3 


1 x 10 4 


B3 


1.0 





30 


3.4 x 10 3 


B4 


0.25 





30 


1.8 x 10 3 



indicated also by a density slope-anisotropy relation proposed 
in 61. 

Concerning the density slope jd, we shall adopt values 
7 d = 1 (NFW profile 01), Id = 1.5 (Moore profile JH) 
and 7d = 0.25 1 for the model B and jd = 12/7 for model A 
OH]. We use NFW or Moore profiles for most of our variants 
of calculation, since these are the most frequently used in the 
literature. 

These density profiles refer to pure dark matter haloes. The 
central part of a halo is likely compressed in the process of 
galaxy formation and baryonic contraction J1,0], which in- 
creases the slope of density profile as well as its normalization. 
To account for adiabatic contraction, we transform f(E, L) 
to /(J, L), where / is radial action - a quantity that is con- 
served during slow change of potential. Hence /(/, L) is con- 
served during compression. As the initial condition for diffu- 
sion equation we take f'(E', L) converted back from /(/, L) 
using now the potential of stars in the bulge and the SMBH. 
(In what follows we omit primes in /' and E'). 

All calculations are performed for the case of Milky Way. 

The normalization of density profile is taken from the con- 
dition that pd(r®) = Pq, where tq = 8 kpc is the distance to 
Galactic center, and pq = 0.3 GeV/cm 3 lIlOll 2 . We summa- 
rize initial model data in Table [I] 

The density profile of bulge stars is taken in the form 

P*(r)=^(fV 7 \ 7,= ( 1 J' r<ri (6) 
\rij [ 2 , r>n 

Here n = 2 pc, i = 6 x 10 4 M Q /pc 3 . Observations sug- 
gest similar broken power-law density profile with the same 
exponents, though with somewhat different values of r\ and 
P*,i ll32ll . The bulge outside SMBH radius of influence is 
therefore taken to be isothermal, with 1-D velocity dispersion 
a = 80 km/s (this is true for inner 10-20 pc, see Fig. 9 in 
14211 ). Relaxation time (O inside ~ 0.2r^ is independent of 



1 This model approximates the inner 10 pc of the Einasto profile in Bertone 
& Merritt IYM . Since it is much more difficult to incorporate profiles with 
variable density slope into our calculations, we have chosen a replacement 
for their "standard halo model" (SHM), that has a very mild cusp with 
average density inside 2 pc equal to that of SHM. 

2 Except for model B4, which is normalized to match SHM of fl3ll . with 
density 40 Mq /pc 3 at 1 pc before contraction. 
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radius and equals 2.6 x 10 9 yr. The velocity distribution of 
stars is assumed to be isotropic. 



B. Diffusion equation and its coefficents 

Gravitational scattering of DM particles by bulge stars leads 
to their diffusion in the {E, L} phase space, which can be 
described by orbit-averaged Fokker-Planck equation 1 35 ] : 



of influence are given below: 



df 

at 



a 

au 



D, 



a/3 



0f_ 

Ota 



Sannlf], (7) 



where £ Q are phase-space variables (usually energy E and an- 
gular momentum L, given per unit mass), Q is the Jacobian, 
D a and D a p are drift and diffusion coefficients, respectively. 
(Given that DM particles have negligible mass, they are not 
subjected to dynamic friction and hence D a = 0). S ann 
is additional loss term due to particle annihilation. We fol- 
low Cohn & Kulsrud lf3~lTl in changing L for dimensionless 
R = L 2 /L 2 , and also change E for dimensionless variable Q 
defined so that Q = 1 over the whole region of interest (both 
the bulge, where E > 0, and the SMBH region of influence, 
where E < 0); Q — ► corresponds to particles very close to 
the SMBH with E — * — oo, Q ~ 1 separates the two regions. 
It can be easily shown that the quantity 



L 2 C (E) I rfi (E) 
3(^)3 



(8) 



satisfies the requirement (here I r> o (E) = I r (E, 0) is the radial 
action of orbit with L = and given E. It can be shown 
that l r {E,R) S3 I r ,o(E) (1 - V~R) HE]). The asymptotic 
expressions for Q are the following: 



Q = - (-E/2a 2 ) 



-3/2 



Q 



3e-y/7T 



exp (3E/2 



E < , (9a) 
E > . (9b) 
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F c is the star distribution function, which is independent of 
energy for 7* = 3/2 in Coulomb potential; In A S3 15 is the 
Coulomb logarithm. (These coefficients may also be derived 
using expressions in IB ill for 7* = 3/2). 

The interesting feature is that Dqq increases with decreas- 
ing R, which means that a particle with small angular mo- 
mentum acquires energy faster (as its orbital pericenter lies 
in the higher density region). Usually this dependence of en- 
ergy diffusion coefficient on angular momentum is ignored, 
and an averaged value is used, which may introduce system- 
atic errors. One of our goals is to check the viability of such 
averaging. 

The expressions for D a p^ for the isothermal bulge are 
somewhat more complicated, so we do not list them here. 
(They can be found in J9[]). The features of expressions are 
similar to ( TTTT i. except that they are multiplied by Q~ 2 / 3 . 

Actually in the calculation we used the coefficients com- 
puted for the exact stellar distribution function, which in turn 
was derived from given density and potential using Eddington 
inversion formula 113511 . These expressions are essentially in- 
terpolating between two limiting cases (actually, close to the 
minimum of two values {D a p c and D a p t b) f° r eacn coeffi- 
cient); the transition occurs at Q ~ 0.05, which corresponds 
to spatial radius of roughly 0.2r^. 



C. Boundary conditions 



From the definition of Q it follows that Q is invariant in 
spherically-symmetric (L = const) and adiabatic (I = const) 
evolution, so that the distribution function in terms of {Q, R} 
is conserved in the process of baryonic contraction and for- 
mation of the black hole. The initial power-law density profile 
corresponds to the following form of the distribution function: 

HQ) = /(JJ)Q-( 6 -T-)/a(4-7 d ) , (io) 

where the dependence on R is the same as in Eq.(l5al>. 

The diffusion coefficients were derived in analytical form 
for the two limiting cases: the bulge (r > r^, or E > 0), 
and the SMBH region of influence, which differ in gravita- 
tional potential and stellar distribution function (details can be 
found in 19J]). The expressions for D a p in the SMBH region 



The diffusion coefficients tend to zero at Q = and R = 1, 
so we do not need any boundary conditions there except regu- 
larity. The presence of black hole induces boundary condition 
at R = R g (Q), where R g = L 2 g /L 2 c , L g = AGM bh /c is the 
minimal angular momentum of a periodic orbit with \E\ <C 
c 2 . The region R < R g is denoted as loss cone. We fol- 
low Lightman & Shapiro [30] in their treatment of the bound- 
ary condition, namely: there exist two absorption regimes - 
empty loss cone, or random-walk capture (for Q < Q cr ), and 
full loss cone, or "pinhole" limit. The key difference between 
them is the ratio of AR, root mean square change of R dur- 
ing one orbital period, to R g , the absorption boundary. In the 
former case, AR <C R g ; particles move in i?,-space by small 
steps in random direction (hence the name), and those getting 
into region R < R g are eliminated in one orbital time. To 
create a flux of particles into the black hole, it is therefore 
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necessary to have gradient in distribution function, df/dR. 
In the opposite case, however, a particle that had R < R g 
at some point of trajectory, can easily diffuse out of the loss 
cone before it actually reaches its orbit pericenter and can be 
removed, and vice versa, particle that initially had R R g 
(but R < Ai?) can be captured by the end of its orbital period. 
So the capture rate depends on amount of particles with small 
R, i.e. on value of / near rather than on its gradient. 
The boundary condition may be expressed as 



= 



R — R, 

where a it taken from Cohn & Kulsrud 1 3 1 



a 



q 



I 0.824^+0.186 9 , q< 1 ; 

\ i , q > i ; 

R.^0 Torb/ Rg 



(12) 

(13a) 
(13b) 



The critical value of Q separating the two regimes, for 
which g = 1, is found to be within SMBH region of influ- 
ence for present-day Milky Way parameters. As seen from 
(1121 1. in the pinhole regime the gradient of / is smaller and 
hence the flux into SMBH is lower [43]. 



IV. SOLUTION OF THE DIFFUSION EQUATION 
A. One-dimensional analysis 

We begin our consideration of DM diffusion with analysis 
of one-dimensional limiting cases. A common approach is to 
assume a logarithmic dependence of / on L (or R), which 
follows from steady-state solution of full 2-D equation lf30L 
[3lll . and solve one-dimensional diffusion equation for energy, 
that includes loss term derived from the solution for angular 
momentum. Actually, this is valid only for quasi-stationary 
solution for stars around black hole, when the loss of stars in 
the BH is balanced by their influx from higher energies. 

In the case of dark matter, the situation is different, since 
there exist no steady state due to the fact that DM particles are 
heated by stars and gain energy rather than lose, and do not 
provide the required influx. So an accurate time-dependent 
solution of full 2-D equation is needed, at least in order to 
verify the validity of reduction to one-dimensional diffusion 
for energy. 

We first discuss the effect of one-dimensional diffusion for 
angular momentum in absence of diffusion for energy. From 
( fTTt ) it follows that Drr w D R for low R. The stationary 

solution of equation = |£ = D-§r [R§r) nas tne f° rm 



f(R) = fg 1 



hi 



R 

Rn 



(14) 



where a is defined in (1 1 3 ab . So for low R the function has 
logarithmic form, as mentioned earlier, and the flux S through 
absorption boundary S = Df g /a is balanced by flux from 
higher R. 



In the model B, f(R) oc R~ fJ at t = 0, and for t > can 
be approximated as 



/,(*)(l + £ln#) , R g <R<Ri(t) 



R-P 



Rx<R^l 



(15) 

So for R < Ri a logarithmic profile has been established, and 
outside i?i the initial value is still conserved; the two profiles 
coincide at R = R%. The values of i?i and f g change in 

time so that the mass loss rate ^ J R f(R, t) dR equals the 
flux to black hole S. This regime holds until R\ reaches its 
maximum of unity (at time t\ ~ 1/D); then the logarithmic 
profile is established for all R, and its amplitude exponentially 
decreases with decay time t — 1/D x (ln(l/i? s ) + a). 

The above approximate arguments give correct qualitative 
results for one-dimensional diffusion along angular momen- 
tum. Global logarithmic profile is established at t > 1 /D and 
decreases in amplitude with characteristic time r (which is 
larger for full loss cone, when a ^S> 1). The loss rate is nearly 
constant (actually, logarithmically depends on f) for t < 1/D 
and decreases exponentially for larger t, proportionally to the 
amplitude of logarithmic profile. 

We now turn to one-dimensional diffusion for energy in the 
absence of loss terms. First we consider the evolution of / in 
the bulge outside SMBH domain of influence. In this region 
Dqq oc Q 4 / 3 , the initial condition has the power-law form 
(fTOl i: f(Q) oc Q~ n ( 5 < n < |), and the one-dimensional 
equation for diffusion along energy can be solved analytically. 
The solution outside SMBH domain of influence may be ap- 
proximated as 



f(Q,t) 



const x t 3 / 2 ™ 

fo(Q) 



Q < Qi(t) , 
Q > Qi(t) . 



(16) 



Qi(t) is defined to match these two asymptotes. This means 
that diffusion effectively blurs / for low Q, driving it towards 
plateau with ever-diminishing amplitude (as a power-law of 
t). 

Then we consider the vicinity of SMBH, where Dqq oc 
Q 2 . The equation then admits solution in the form f(Q) oc 
Q~ n exp(— t/r) + const, where r is the relaxation time (as 
said above, it is independent of energy in the case considered, 
but still depends on R: relaxation goes on faster for smaller 
angular momenta). If we neglect the dependence of t on R, 
we may come to the following conclusion: as long as t < t, 
the evolution of / (and hence density) is only in amplitude 
which drops exponentially, not the shape. When the exponent 
becomes significant (or t > r), the function tends to constant 
value which is set by the outer boundary condition (where it 
is determined by the solution outside SMBH domain of influ- 
ence described above). Since the relaxation time in the inner 
area is always shorter than in bulge, this constant (with respect 
to Q) value is determined solely by the diffusion in bulge and 
is decreasing in time. The constant value of / implies the 
density profile p oc r~ 3 / 2 , as mentioned in l37ll . (This is dif- 
ferent in the case of stars, where the effect of dynamic friction 
should be included in the equation, and the corresponding so- 
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TABLE II: Variants of calculation. Initial halo models are described in Table U Mi„u is the initial DM mass within rn = 2 pc, M cap t - 
amount of DM captured by SMBH (in parentheses - from r < r^), M evap - DM mass evaporated from r < r^, M ann - annihilated mass (all 
of these - for present time, t — 10 10 yr, in Mq). The final two columns give logarithm of astrophysical factor J in annihilation flux, averaged 
over solid angle 10 -3 sr (and 10 -5 sr in parentheses), for t — and t — tnubbie, respectively. 

All variants except 7 have Mbh = const = 3 x 10 6 Mq, in var.7 Mbh rises from 3 x 10 4 M© to its present-day value as \/t. All variants 
except 6 are solved using two-dimensional equation; var.6 uses simplified one-dimensional diffusion equation for energy as described in 
Section lTV Al In variants 2 and 5 maximal annihilation cross-section {u a v x } = 3x 10~ 26 cm 3 /s and particle mass m x = 50 GeV was used, 
in others annihilation was not included. 



Var model, features 


Mi n it 


M capt {Mt pt ) 






log X(Js),*=0 


-"- t — tHub 


1 


B 1 Moore, reference var. 


1.32 x 10 6 


1.78 x 10 5 


(1.00 x 10 s ) 


5.89 x 10 5 




15.61 (17.61) 


6.20 (8.19) 


2 


B 1 maximal annihilation 


1.32 x 10 6 


1.64 x 10 5 


(8.69 x 10 4 ) 


5.77 x 10 5 


3.04 x 10 4 


15.61 (17.61) 


6.13(8.12) 


3 


B2 radial anisotropy 


9.79 x 10 5 


1.94 x 10 s 


(8.55 x 10 4 ) 


3.98 x 10 s 




14.54 (16.54) 


5.93 (7.93) 


4 


B3 NFW, reference var. 


3.05 x 10 5 


3.86 x 10 4 


(2.01 x 10 4 ) 


1.28 x 10 s 




13.56 (15.56) 


5.05 (7.04) 


5 


B3 maximal annihilation 


3.05 x 10 5 


3.73 x 10 4 


(1.88 x 10 4 ) 


1.27 x 10 5 


2.72 x 10 3 


13.56 (15.56) 


5.00 (6.99) 


6 


B3 1-D diffusion 


3.05 x 10 5 


4.05 x 10 4 


(2.08 x 10 4 ) 


1.31 x 10 5 




13.56 (15.56) 


5.03 (7.02) 


7 


B3 M bh oc s/i 


2.69 x 10 5 


2.17 x 10 4 


(8.85 x 10 3 ) 


1.36 x 10 5 




14.22 (16.22) 


5.02 (7.01) 
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B4 Einasto profile 


1.47 x 10 5 


1.77 x 10 4 


(8.49 x 10 3 ) 


5.74 x 10 4 




12.13 (14.13) 


4.56 (6.54) 
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Al Gurevich-Zybin profile 


2.17 x 10 6 


6.75 x 10 5 


(2.68 x 10 5 ) 


7.22 x 10 s 




12.26 (14.26) 


6.56 (8.56) 



lution has the Bahcall-Wolf density profile p oc r~ 7 / 4 ). How- 
ever, the profile is not stationary and its amplitude decreases 
in time (as noted in Q3S0 ), which is attributed to ongoing heat- 
ing of dark matter particles by stars, continued mainly outside 
rh after one relaxation time. The amplitude then decreases as 
power-law of time: p ex £- 3 / 2 ™ = t-(6-7)/(8-2 7 )_ 



We want to draw attention to the importance of considering 
both regions (bulge and SMBH sphere of influence) in order 
to obtain the density profile and annihilation signal evolution 
for t ~ t r center . For both cases t ^> t r and t -C i r analytical 
estimates can be obtained, but for transition period t r < t < 
4t r exact calculation is required. Specifically, Galactic center 
is just in this category. 



Thus far we have considered only one-dimensional approx- 
imations, with the diffusion along angular momentum caused 
by capture of DM particles by the black hole, and diffusion 
along energy caused by heating on stars. A common approach 
to solution of the two-dimensional equation is to consider one- 
dimensional diffusion along energy with additional loss term 
due to capture of particles by black hole, caused by diffusion 
along momentum. The loss term is taken from the stationary 
solution of diffusion equation for momentum, and introduced 
into the equation for energy as an additional term S. The equa- 
tion itself is written for averaged f(Q) = J Q f(Q,R)dR. 
The magnitude of the loss term S is taken from the boundary 
condition at the black hole: S = D f g /a, where f g is related 
to / by O. The diffusion coefficient Dee also should be 
averaged over R, despite the fact that it significantly varies 
with R (see (fTTb)). The validity of this reduction should 
therefore be verified by comparison with solution of full two- 
dimensional equation, which was one of aims of the present 
study. 



B. Two-dimensional equation: variants of calculation 

In order to address different issues mentioned in SectionlTll 
we have performed several variants of calculation, which dif- 
fer in initial conditions and in account for distinct factors of 
evolution. All of them are listed in Table iHl 

The full two-dimensional diffusion equation (0 was inte- 
grated on rectangular grid with variable space- and time-steps. 
The boundary condition at R = R g was treated in the same 
way as in IBlll (section Vc), i.e. for R g less than ARi (the size 
of first cell), flux S through absorption boundary was consid- 
ered as flux through R = 0, which was related to the value 
of / in the center of the cell via (O with S = Df g /a. The 
spatial density profile was recalculated every few timesteps 
and used to compute orbit-averaged annihilation rate at each 
{Q, R} in variants 2 and 5. 

The variants are chosen to investigate the relative role of 
different factors. Five different initial models were chosen 
(Sec. IHL AD . The model B3 (NFW profile with isotropic ve- 
locity distribution) is a basic model to study the influence the 
one-dimensional approximation (Var.6) and growing BH mass 
(Var.7), which are to be compared with reference Var.4. From 
the existence of quasars at z w 6 B44I1 we know that SMBHs 
were already massive at very early times, but the situation is 
also possible when the black hole grows slowly from small 
initial seed due to accretion of stars and gas. The growth law 
is taken to be Mhh oc \fi [45]; the SMBH radius of influence 
and the total gravitational potential are changed accordingly. 
In other variants Mhh was held constant, and gravitational po- 
tential remained fixed. 

The influence of velocity anisotropy is studied in Var.3 for 
Moore profile (from the (3 — 7 relation of Hansen & Moore 
114111 it is naturally to expect more radial anisotropy for steeper 
profiles), which should be compared to Var. 1 . 

Additionally, annihilation was accounted for in Vars.2 and 
5, with maximal annihilation cross-section of (a a v x ) = 3 x 
10~ 26 cm 3 /s and particle mass m x = 50 GeV, same as in 
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111 Oil . These should be compared to Vars. 1 and 4, respectively. 
In other variants annihilation was disregarded, but predictions 
of 7-ray flux were made (see below). 

The prospects of indirect DM detection are related to regis- 
tration of 7-rays from the DM annihilation in the most dense 
regions. One of such sites is the Galactic center, where in- 
deed a 7-ray source was detected by several ground-based 
Cherenkov telescopes such as H.E.S.S. S, VERITAS 63] , 
CANGAROO B and MAGIC & as well as space tele- 
scope EGRET |5(J]. 

The 7-ray flux is usually split into two factors, one of them 
depending on particle physics parameters and the other (J - 
astrophysical factor) is determined by DM density distribution 
& 



$ ~ $ 



>aup/ 



3 x 10- 26 cm 3 /s 



lOGeV 



(17) 



where <& = 5.6 x 10 8 cm 2 s 1 , and Jao is dimensionless 
integral of squared density, averaged over a solid angle AS1: 



1 

Ml 



dl 



/9 2 (r)27n/; dtp 



(18) 



K~ x = (8 kpc)(0.3 GeV/cm 3 ) 2 , and integral is taken over 
line of sight and the solid angle corresponding to telescope 
resolution. (For simplicity, the point-spread function is taken 
as Heaviside step function). 

We calculate values of J5 and J3 corresponding to Afi = 
10~ 5 and 10~ 3 . The former is the approximate resolution 
of GLAST lf52tl and modern Cherenkov telescopes such as 
HESS, the latter is the resolution of EGRET [50]. The spa- 
tial resolution is 15 and 150 pc, correspondingly. It is obvious 
that for a power-law density profile p oc r~ 7 with 7 > 1.5 the 
main contribution in ( fT~8T > comes from the inner boundary, and 
for 7 < 1.5 the situation is inverse. Since we have adiabat- 
ically compressed dark matter halo, its density slope 7' after 
contraction is always greater than 1.5 (7' = (6 — 7)/(4 — 7), 
where 7 is initial slope). Therefore, the flux is determined by 
dark matter distribution in the vicinity of black hole, and the 
additional contribution from the dark matter in bulge (r 3> r^) 
is negligible. This is in contrast with the conclusion of l39ll . 
where the authors did not use adiabatically contracted dark 
matter profile and hence the contribution from bulge was sig- 
nificant because 7 < 1.5. 

Another observational constraint on dark matter distribu- 
tion may come from the study of stellar orbits which measure 
the enclosed mass and may help to determine the fraction of 
dark mass near the SMBH JH]. 

The results of calculation are presented in TableM The five 
initial models differ in normalization, i.e. in the initial dark 
matter mass Mj n # within r^. Since the dynamics of dark 
matter is governed by linear equations (except for annihila- 
tion), the results may well be scaled to match each other. The 
similarity between different models then becomes apparent. 
Approximately 10-15% of Mi n n is captured by black hole by 
t = 10 10 yr (M capt ); half of this mass comes from the SMBH 
sphere of influence (r < r^, M^ apt ), About 40% of Mi n n has 
been evaporated from the central region. 
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FIG. 1: Density profiles for variants 4 (without annihilation, solid), 
5 (with maximal annihilation, dot-dashed), and for pure annihila- 
tion without diffusion (dotted lines). Top curve in each series is 
for t = 2.5 Gyr, approximately one relaxation time; bottom is for 
t — 10 Gyr. Initial adiabatically compressed profile is shown by 
long-dashed line. 

Without diffusion, the density profile would be a broken power-law 
with p oc r~ 0,5 for r < r a , and p oc r -2,33 for r > r a (initial pro- 
file); r a is the annihilation radius at which the initial density equals 
p a i© which is set at 1.4 x 10 8 M & pc -3 by now. 
Without annihilation the density would drop several times by one re- 
laxation time in the whole region r < rn (top solid line), and drop 
dramatically by now (bottom solid line); since it becomes shallower 
than par"' 5 (which is sketched as dashed line) inside r^, the 
main contribution to annihilation factor J comes from r ~ r^. 
Combining together the effects of annihilation and diffusion (dot- 
dashed lines), we get profiles that are similar to pure diffusional at 
large radii and are cut by annihilation at very small radii. The J fac- 
tor by now is essentially the same as without annihilation and is not 
affected by this cutoff, but for t = t r the difference in J between 
Var.4 and Var.5 is still large (see Fig©. 



The difference between variants with the same initial model 
are significant only for t < At r — 10 10 yr, which is by coin- 
cidence the total time of evolution. t r ^ is the relaxation 
time inside and is independent of radius. So the current 
state (and further evolution) is basically independent of initial 
conditions, although if t r was, say, twice as large, then the 
difference would be much more prominent by present time. 

First we consider basic variant 4 without annihilation. 
Since the diffusion along energy goes on faster for low an- 
gular momenta, the particles with both low energy and low 
angular momenta are quickly heated up or captured by black 
hole, therefore the density in the immediate vicinity of the 
black hole quickly drops. However, by one relaxation time 
it decreases roughly by the same factor of few in the whole 
region r < 0.27-^ (Fig.Q] upper solid), except for the very in- 
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FIG. 2: Evolution of annihilation factor J 5 for models 4 (with- 
out annihilation, solid), 5 (maximal annihilation, dot-dashed), 7 
(Mth oc \ft, dashed), and for pure annihilation without diffusion 
(dotted lines). 

The latter one is a power-law in time (J oc t~ 5//7 ); the case without 
annihilation has a sharp decline at t > t r ~ 2.5 Gyr which corre- 
sponds to transition from a steeper than r~ 1,5 slope before relaxation 
to a shallower slope when the value of J is dominated by density at 
r ~ r/j and changes very weakly in time since t > 4t r . 
The model 5 with both annihilation and diffusion included is roughly 
an interpolation between these two extremes: before relaxation the 
value of J is dominated by density at small radii (~ r a ), and after it 
is identical to the model 4 without annihilation. The time after which 
values of J for models 4 and 5 start to be equal happens to coincide 
with current time, which is by no means an rigorous result. 
The model 7 has much shorter relaxation time so the sharp decline 
occurs earlier, and since then the evolution is close to self-similar 
with very slow decline of J in time. 



tier radii. By this time the density profile is still steeper than 
r -1,5 inside rh, and hence the integral J p 2 (r)r 2 dr which de- 
termines the J-factor in ( fl"8l depends on the low-r density be- 
havior. Then the heating starts to dominate, and within few t r 
the distribution function and density tend to nearly universal 
form: for r < O.lr/j density p ~ r -15 , near r ~ r h density 
profile is a little bit shallower, and further away it is identical 
to initial profile. So most of J is gained at r ~ ru, and hence 
the difference between variants 4 — 7 almost vanishes. The 
density profile inside then drops nearly self-similarly (if 
we continue integration to 10 11 yrs, density is reduced about 
an order of magnitude everywhere at r < r/, while retaining 
its form). It is important to note that after this regime has 
been established, diffusion outside SMBH sphere of influence 
plays the key role. 

In the models with initial dominance of low angular mo- 
mentum orbits (B2 and Al) the initial flux of dark matter into 
the SMBH is larger, but then as the it'-diffusion establishes 
near-logarithmic profile (at t > t r ), the memory of initial con- 
ditions is erased. 

The inclusion of self-annihilation (in Vars. 2 and 5) affects 
only the inner region in the manner that the density does not 
exceed the "annihilation weak cusp" (not plateau). Conse- 



FIG. 3: Evolution of dark matter mass captured by black hole, for 
models 3 (dotted), 1 (long-dashed), 4 (solid) and 7 (short-dashed 
line), from top to bottom. Models 1 (Moore profile) and 4 (NFW pro- 
file) differ in initial dark matter mass and hence the same difference 
in captured mass, with similar dynamics. Model 3 is Moore profile 
with initial radial anisotropy: initially the captured mass grows faster 
because of greater abundance of particles with low angular momenta; 
however, after one relaxation time the difference in initial conditions 
almost disappears and the capture rate and total mass becomes sim- 
ilar. Model 7 features small initial seed black hole mass and shorter 
relaxation timescale, hence the capture rate is lower and heating goes 
on faster; by now the total captured mass is smaller. 



quently, the J value diminishes faster than without annihi- 
lation (Fig. |2]), but after t > t r the heating starts to prevail 
over annihilation even in the very inner region, and everything 
again goes as in basic variant. 

The variant 7 which features a black hole growing from 
small initial seed differs in the way that relaxation time inside 
rh is less and a nearly self-similar relaxed solution is estab- 
lished earlier. In this case the total dark matter mass captured 
by the black hole is also less than in the case of constant black 
hole mass (Fig. O. In all variants the total captured mass is 
much less than Mbh, so the neglect of the absorption on the 
black hole mass is justified. 

Finally, comparison of one-dimensional solution (Var.6) to 
the two-dimensional one (Var.4) reveals very little difference 
in either final result or dynamics. This is because i?-profile of 
distribution function at fixed Q is reasonably well described 
by a logarithm, failing only at very early (t <C t r ) or interme- 
diate (t ~ tr) times (Fig.|4|i. Therefore, the commonly used 
approach of reduction to one-dimensional diffusion equation 
is a good approximation, at least for the class of initial distri- 
bution functions that have weak dependency on R (the case 
of velocity isotropy - model B3 - being in this class). More- 
over, the evolution of i?-averaged distribution function is qual- 
itatively consistent with simple analytical description of Sec- 
tion |TV_A] (Fig . E) . 
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FIG. 4: Distribution function f(R) for fixed Q: solid lines - Var.4, 
dotted lines - Var.6 (1-D treatment). Curves A and C are for Q = 
10~° (corresponding to circular orbit with r c ; rc = lCP^r^, R g — 
0.0011, and a = 1.6 X 10" 3 (Eq©), curve B is for Q = 10" 3 
(r circ = 0.1r h , R g = 1.1 x 10~ 5 , and a = 0.56). Curves A and C 
are for t — 0.01i r , curve B is for t = t r . 

The difference between solution of exact two-dimensional equation 
and simplified one-dimensional treatment, where profile is exactly 
logarithmic (Eqll4b. is rather moderate: for early times (curves A and 
B) the logarithmic profile in the exact solution is not yet established 
for all R, and for high R the function is close to initial constant value; 
while for intermediate times (t ~ t r , curve B) the difference is due to 
intrinsically two-dimensional character of evolution. (For t = 0.1t r 
or t = 10t r the corresponding curves would be very close). In all 
cases, the behavior of f(R) near loss-cone boundary R g is very close 
to logarithmic; the difference from 1-D case is only in amplitude. 
The curves A and B illustrate random-walk regime, or empty loss 
cone; curve C suggests that in pinhole regime (a 2> 1) f(R) is 
approximately constant. 

V. DISCUSSION AND CONCLUSIONS 

We have considered the evolution of dark matter distribu- 
tion in the Galactic center caused by gravitational scattering 
on stars of galactic bulge, absorption by the supermassive 
black hole and annihilation. The spherically symmetric evo- 
lution is described by two-dimensional Fokker-Planck equa- 
tion in phase space of energy and angular momentum. One 
of main goals was to investigate how the difference in initial 
conditions and in various dynamical factors affects the final 
result. 

It turns out that at least for our Galaxy there has been 
enough time to eliminate most of the differences and to es- 
tablish a nearly universal distribution and density profile with 
the following features: for r < 0.1 w 0.2 pc the density is 
close to power-law p oc r~ 7 with slope j as 1.5 or less, far- 
ther out the slope becomes shallower, and even farther it joins 
the initial profile which is already steeper than r~ 15 because 
of adiabatic compression in the bulge. Therefore, as regard- 
ing the possible annihilation signal, after approximately four 
relaxation times it is determined by DM distribution near Th 
and depends only on initial profile normalization, not the de- 
tails of evolution. This explains the low scatter in J values 
in the last column of Table UH the most probable range for 
J5 being 10 6 — 10 8 . (This is in moderate agreement with re- 
sults of Bertone & Merritt [13], though in our calculations the 
difference between models with no annihilation and maximal 
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FIG. 5: Evolution of i?-averaged distribution function in Var.4 for 
different t (from top to bottom): t = 0, t = 2.5 x 10 9 yr » t r , 
t = 5 x 10 9 , t = 10 10 and t = 10 11 yr. Initially / oc Q" 5/9 ; by one 
relaxation time / drops by an order of magnitude nearly uniformly 
inside SMBH sphere of influence, where t r is constant; by 4t r it is 
almost flat for Q < 1 and continues to drop as t , as described 
in Section |IV A| The slope Q~ 5 ^ 9 corresponds to density profile 
p oc r~ 7 ' 3 inside rn and r _5//3 outside; flat distribution - to r~ 3 ' 2 
and r , correspondingly. The top axis shows r c - radii of circular 
orbits for given Q, in units of ru- 

annihilation is much smaller). The low scatter in J means that 
constraints on the microphysical nature of dark matter follow- 
ing from Eq.dT7t will be more robust, if the 7-ray flux detected 
from the Galactic center is at least partially attributed to dark 
matter annihilation 11521 15411 . 

The different initial models also demonstrate similar pro- 
portions of dark matter mass that has been captured by the 
black hole (~ 10 — 15%) and evaporated from the central 
2 pc (~ 40% of the initial DM mass within 2 pc). A num- 
ber of previous studies suggested that the DM may constitute 
a significant fraction of the total SMBH mass. However, in 
these papers either the loss cone was assumed to be always 
full II551I56I1 . or the diffusion along energy axis (heating) was 
disregarded it @], which led to higher estimate of the cap- 
tured mass compared to the present study. It remains in ques- 
tion whether a non-spherically-symmetric initial distribution 
could sustain a full "loss cone" due to chaotic randomization 
of orbits l55l H3]. This issue will be addressed in a future 
paper. 

We found that the one-dimensional diffusion approximation 
works good for models with / initially independent on angular 
momentum. 

Comparison with paper 13811 where 1-D equation was 
solved shows qualitatively similar behavior of i?-averaged 
distribution function and density evolution: for t <t r the den- 
sity drops self-similarly in the SMBH region of influence, and 
for t t r it tends to r~ 3 / 2 form. The quantitative difference 
might be attributed to different stellar density profile outside 
~ O.Itv the evolution in the transition spatial (r ~ 77J and 
temporal (t ~ t r ) regions essentially depends on the magni- 
tude of diffusion coefficients which are determined by stellar 
density. 

There are other targets of indirect dark matter search, ex- 
cept the Galactic center, such as nearby galaxies or globular 
clusters M58I1 . The importance of the processes discussed in 
this paper on the 7-ray flux from these other targets depends 
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on whether they are dynamically old, i.e. is their relaxation 
time in the very center less than Hubble time. For exam- 
ple, dark matter dominated dwarf spheroidal galaxies such 
as Draco have relaxation times well exceeding 10 10 yr, so 
the scattering off by stars is unimportant, and the flux is de- 
termined by initial conditions and annihilation cross-section 
ll59tl . On the other hand, globular clusters such as Gl in An- 
dromeda Galaxy are dynamically old, and are expected to be 
on the relaxed stage of evolution and therefore have less scat- 
ter in J, hence they are also promising in constraining dark 
matter parameters l60ll . (However, in this case one should 
account for evolution of stellar distribution as well). Finally, 
galaxies on the intermediate stage of evolution such as M32 



which has central relaxation time of order 2 — 3 x 10 9 yr l6lll 
may exhibit very different values of J (a drop by many orders 
of magnitude during several t r is seen on Fig. [2j. 

The concentration of dark matter in the center of our Galaxy 
is a promising target of indirect dark matter detection, and 
its dynamical evolution is crucial for constraining dark matter 
parameters. 
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